function [Trc,Tspc,u,u_10percent,c_hansen,c_hansen_10percent,naive_p,naive_10percent,best,besthansen]=bsds_ralf_ess(bench,models,B,w,boot,number)
% PURPOSE:
% Calculate Whites and Hansens p-vals for outperformance
% 
% USAGE:
% [c,u,l]=bsds(bench,models,B,w,boot)
% 
% INPUTS:
% bench  - The return series form the benchmark model
% models - The return series from each of the models used for comparrison
% B      - Number of Bootstrap replications
% w      - Desired block length
% boot   - 'STATIONARY' or 'BLOCK'.  Stationary will be used as default.
% 
% OUTPUTS:
% c - Consistent P-val(Hansen)
% u - Upper P-val(White)(Original RC P-vals)
% l - Lower P-val(Hansen)
% 
% COMMENTS:
% 
% 
% Author: Kevin Sheppard
% kksheppard@ucsd.edu
% Revision: 2    Date: 12/31/2001

% Revision: a few times by Ralf Steinhauser in 2006


[t,k]=size(models);


if size(bench,2)~=1
    error('Only 1 benchmark allowed');
end


if t~=size(models,1)               
    error('Data and Models must have the same length');
end

if nargin==6
    if strcmp(boot,'STATIONARY')                
        p=1/w;
        [bsdata]=stationary_bootstrap([1:t]',w,B); 
    else
        [bsdata]=block_bootstrap([1:t]',w,B);
    end
end

if nargin<6
    [bsdata]=block_bootstrap([1:t]',w,B);
end


%OK now we have the bootstraps, what to do with them?
diffs=(repmat(bench,1,k)-models);       

%The the consistent
%need to bootstrap the data here as well
A=zeros(k,1);
 
for i=1:k
    %[bsdata2]=stationary_bootstrap([1:t]',1/w,B);         
    workdiffs=diffs(:,i);
    temp=workdiffs(bsdata);                 
    temp=mean(temp);
    A(i)=1/4*t^(0.25)*sqrt(sum((t^(0.5)*temp-t^(0.5)*mean(diffs(:,i))).^2)/B);  
    std(i)=sqrt(sum((t^(0.5)*temp-t^(0.5)*mean(diffs(:,i))).^2)/B);
end

stat=max(mean(diffs));
best=find(mean(diffs)==max(mean(diffs)));

stathansen=max(mean(diffs)./std);
besthansen=find((mean(diffs)./std)==max(mean(diffs)./std));


g=(mean(diffs)>=-A').*mean(diffs);

%We now need to build up our matrix of means for each series and each bootstrap, k x B
perf=zeros(B,k);
for i=1:k
    workdata=diffs(:,i);
    perf(:,i)=mean(workdata(bsdata)-g(i))';
end
perfc=max(perf,[],2);              
c=mean(perfc>stat);

perfch=max((perf./repmat(std,B,1)),[],2);
c_hansen=mean(perfch>stathansen);

if nargout>1
    %Then the upper
    g=mean(diffs);
    perf=zeros(B,k);
    for i=1:k
        workdata=diffs(:,i);
        perf(:,i)=mean(workdata(bsdata)-g(i))'; 
    end
  
    perfu=max(perf,[],2);
    u=mean(perfu>stat);
 
    perfo=max((perf./repmat(std,B,1)),[],2);
    u_hansen=mean(perfo>stathansen);
end

if nargout>2
    % First the lower
    g=max(0,mean(diffs));
    perf=zeros(B,k);
    for i=1:k
        workdata=diffs(:,i);
        perf(:,i)=mean(workdata(bsdata)-g(i))';
    end
    perfl=max(perf,[],2);
    l=mean(perfl>stat);
    
    perfo=max((perf./repmat(std,B,1)),[],2);
    l_hansen=mean(perfo>stathansen);
end

c_distribution=sort(perfc);
c_10percent=c_distribution(0.9*B)*t^(0.5);

c_distribution=sort(perfch);
c_hansen_10percent=c_distribution(0.9*B)*t^(0.5);

u_distribution=sort(perfu);
u_10percent=u_distribution(0.9*B)*t^(0.5);

l_distribution=sort(perfl);
l_10percent=l_distribution(0.9*B)*t^(0.5);

Trc=stat*t^(0.5);
Tspc=stathansen*t^(0.5);

number=find(mean(diffs)==max(mean(diffs)));
workdata=diffs(:,number);
naive_distribution=sort(mean(workdata(bsdata)-mean(diffs(:,number))));
naive_10percent=naive_distribution(450);


naive_p=mean(naive_distribution>stat);

